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Robust oscillations in SIS epidemics on adaptive networks: 
Coarse-graining by automated moment closure 
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Abstract. - We investigate the dynamics of an epidemiological susceptible-infected-susceptible 
(SIS) model on an adaptive network. This model combines epidemic spreading (dynamics on the 
network) with rewiring of network connections (topological evolution of the network). We propose 
and implement a computational approach that enables us to study the dynamics of the network 
directly on an emergent, coarse-grained level. The approach sidesteps the derivation of closed low- 
dimensional approximations. Our investigations reveal that global coupling, which enters through 
the awareness of the population to the disease, can result in robust large-amplitude oscillations of 
the state and topology of the network. 
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Introduction. — Over recent years the physics of net- 
works has received much attention [1-3]. Research has 
focused on the one hand on dynamics of networks, i.e., 
network growth and rewiring [4,5], and on the other hand 
on dynamics on networks; the coupling of individual dy- 
namical systems according to a static network topology 
[6,7]. Only recently have these two distinct strands of 
network research been brought together in the study of 
adaptive networks, which combine topological evolution of 
the network with dynamics on the network nodes [8-14]. 
Adaptive networks have been shown to exhibit a num- 
ber of new phenomena, including robust self-organization 
towards dynamical criticality [8], the formation of com- 
plex global topologies based on simple local rules [13,15], 
the emergence of new bifurcations and phase transitions 
involving local as well as topological degrees of freedom 
[12,13], and finally a spontaneous 'division of labor' in 
which an initially homogeneous population of nodes self- 
organizes into functionally distinct classes [10] . 

A prominent application of networks is epidemiology. 
In an ever more populated and tighter connected world 
epidemic diseases are once again becoming a major threat 
[16]. In the past decades we have witnessed the return of 



old diseases [17] as well as the emergence of new ones [18] 
at an accelerating pace. Conceptual mathematical models 
can provide fundamental insights on the underlying mech- 
anisms that govern epidemic dynamics. 

One of the most simple conceptual models of epidemics 
is the susceptible-infected-susceptible (SIS) model [19], de- 
scribed in more detail below. If considered on a static 
network, varying the parameters in this model reveals at 
most one dynamical transition. This transition is continu- 
ous and corresponds to the epidemic threshold: the point 
in parameter space beyond which the disease can invade 
the network. If one enables the individuals in the model 
to avoid contact with infected ones by altering their lo- 
cal interaction topology the system turns into an adaptive 
network. It has recently been shown that this natural ex- 
tension of the spatial SIS model gives rise to more complex 
dynamics, including bistability, hysteresis and a narrow 
region of oscillatory dynamics [13]. 

In order to understand the subtle interplay between 
local dynamics and topological evolution that character- 
izes adaptive networks the application of tools from dy- 
namical systems theory is desirable. While these meth- 
ods arc in general only applicable to relatively low- 
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dimensional equation-based models, previous results [13, 
20] suggest that often a small number of topological de- 
grees of freedom suffices to characterize the dynamics of 
the individual-based model on an emergent level. This 
means that a small number of observables can describe 
the macroscopic (system-level) state of the network and 
thus a low-dimensional description at this level is feasible. 
In [13] a low-dimensional set of system-level equations of 
motion for the adaptive SIS model was derived by means 
of a moment expansion and subsequent moment closure 
approximation. While this approach has yielded satisfac- 
tory results, it involves a strong homogeneity assumption 
and therefore can fail in certain topologies. 

In this Letter we extend the previous work in two ways. 
First, we consider the effect of a variable level of awareness 
to the disease. This awareness is shared instantaneously 
among the population and therefore acts as a global cou- 
pling. Second, we propose and implement an alterna- 
tive computational approach, based on the coarse-grained 
'equation-free' modeling and analysis framework [21]. In 
this approach the analytical moment-closure approxima- 
tion is replaced by a numerical procedure which automat- 
ically extracts the effects of appropriate closure. This ap- 
proach enables us to investigate the emergent dynamics of 
complex networks with well-established tools of dynamical 
systems theory without explicitly deriving a closed analyti- 
cal description on this level. Our investigation reveals that 
the global coupling strongly extends the parameter region 
in which oscillatory dynamics can be observed. It thereby 
gives rise to persistent large-amplitude oscillations affect- 
ing the prevalence of the disease as well as the topology of 
the network. 

The adaptive SIS model. — We consider a network 
with a fixed number of nodes N and undirected links L. 
Every node represents an individual, while links represent 
social contacts between individuals. An individual can be 
either infected (I) or susceptible (S). We denote the links 
between individuals as SS-links, Si-links or II-links accord- 
ing to the states of the nodes that they connect. In every 
time step infected individuals recover with a probability r, 
becoming susceptible. For every Si-Link there is a proba- 
bility p that the susceptible individual becomes infected. 
In addition to these standard rules of SIS models, we allow 
susceptible individuals to avoid contact with infected in- 
dividuals by altering their local interaction topology: For 
every Si-Link there is a probability w that a rewiring event 
occurs in a given time step. In a rewiring event, the suscep- 
tible individual cuts its link to the infected individual and 
establishes a new link to another susceptible individual. In 
contrast to [13] we do not fix the rate w, but assume that 
it changes according to the individual's awareness of the 
disease. In human populations with access to mass media 
it is reasonable to assume that information on the dis- 
ease is instantaneously disseminated across the network. 
We therefore set w = wop, where p = i/N , is the infected 
fraction of the population and wq is a constant parameter. 



Coarse graining. — It has been shown in [13] that 
the system-level dynamics of sufficiently large networks 
(here N = 10 5 , L = 10 6 ) can be captured by the number 
of SS-links Zgs, the number of II-links hi and the number 
of infected individuals i. The number of susceptibles s and 
the number of Si-links Zsi then follow from the conserva- 
tion laws s + i = N and Zss + Isi + hi = L. A moment 
expansion reveals that the dynamics of the collective state 
variables s, Iss, hi can be described by 

■±i = p l sl - ri 

ihi = Pisi(^ + i)-nm (1) 

&lss = ri l -f-pl SI l -f^+wl SI 

where Zisi and Zssi denote the number of I-S-I and S-S-I 
chains in the network, respectively. 

Note that we have so far only assumed that the net- 
work is large, so that the state can be described by con- 
tinuous variables and the effect of individual-based (de- 
mographic) stochasticity can be neglected. In particular 
no strong homogeneity has been used since many effects 
of non-homogeneous topology are captured by taking the 
densities of links lss,hn,hi and chains Zisi,Zssi into ac- 
count. However, Eq. (1) docs not yet constitute a closed 
model, as we haven't specified how Zisi and Issi are deter- 
mined. 

In Rcf. [13] the moment-closure approximation Zisi = 
^si/ s '^ssi = 'ss'si/2s was used to close the model. This 
approximation introduces a homogeneity assumption as 
it requires that the SI- and SS-links are homogeneously 
distributed among the possible three node chains. While 
the moment closure approximation has yielded good re- 
sults in [13] it has to fail in networks with slowly decaying 
degree distribution, such as scale free graphs. 

In this Letter we avoid the moment closure approxima- 
tion: We do not attempt to derive analytical expressions 
for Zisi and Zssi m order to close Eq. (1). Instead, we 
specify a numerical procedure that generates the effect of 
appropriate closure terms on-demand from short bursts of 
numerical simulation. Therefore the major technical chal- 
lenge in this Letter is to determine the correct value of Zisi 
and Zssi for a given set of state variables i, Iss, hi- 

In general one can image a large number of candidate 
network topologies that agree with the desired values of 
i, his, hi- Thus, for a given set of i, Iss, hi, many different 
values of Zisi and Zssi are in principle possible. However, 
most of the candidate topologies will not arise in the long 
term dynamics of the network. In fact one finds that the 
Zisi and Zssi that are consistent with given i,Zss,Zn are 
restricted to a small set in the long-term dynamics. 

Consider the following: Large adaptive networks can 
easily have millions of topological degrees of freedom. The 
assumption that the state of such a network can be char- 
acterized by a low number of variables implies in general 
that there is a time scale separation between a few slow 
variables and the remaining fast degrees of freedom. The 
system will then approach a slow manifold, on which the 
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Second Moment Li 

Fig. 1: (Color online) Schematic of the lifting procedure. A 
network with initially arbitrary third moments is settled to 
the slow manifold (dotted line) on which the consistent val- 
ues encountered in the long-term dynamics are located. For 
this purpose we run a short simulation (solid line) and then 
reset the first and second moments while retaining the third 
moments (dashed line). Inset: Illustration of coarse projec- 
tive integration. The necessary information for long projective 
leaps (dashed) is extracted from short bursts of individual- 
based simulation (solid). 

long term dynamics takes place. The slow manifold can be 
parameterized entirely by the slow variables, as all other 
topological degrees of freedom are slaved to the slow vari- 
ables. Thus a set of slow variables determines the state of 
the system. This implies that all other degrees of freedom 
can be written as a function of the slow variables. 

In the system considered here the time scale separation 
between the slow variables i, l$s, hi and fast variables, such 
as Issi and Zisi, is finite and thus the reasoning given above 
is only approximately true. Nevertheless, the values of Zisi 
and Zssi that can possibly appear in the long term dynam- 
ics (while a given set of i,lss, hi is observed) are restricted 
to a relatively narrow band. This band is simply a noisy 
graph of the slow manifold. If the time scale separation is 
sufficiently large, the band of possible Zisi or Zssi is nar- 
row enough to treat it effectively as a single value. In 
the following we denote such a value of Zisi or Zssi that 
is consistent with a given set of i, ZsS; hi i n the long term 
dynamics simply as a consistent value. 

Let us now describe an algorithm that computes the 
consistent values of Zisi and Zssi for a given set of state 
variable values i*, Zg S , Zjj. In order to account for the com- 
binatorial effect of the slow variables, it is advantageous to 
work with the normalized third moments cisi = sZisi/Zsi 2 , 
cssi = hsis / (hshi) ■ We start by randomly generating a 
network at the desired operating point i* , Zjj, Zg g . The 
third moments Zisi and Zssi i n a random graph will not, 
in general, be consistent, that is, they are not on the slow 
manifold of the system. In order to "settle" the network 
to the slow manifold, where the consistent values are lo- 
cated, we proceed as follows: We run a short simulation 
of j steps. Over the course of the simulation the sys- 
tem approaches the (attracting) slow manifold; the val- 
ues of the third moments at the end of the simulation 




Fig. 2: A slice of the slow manifold (i = 97693, l BB = 884707). 
Shown are the closure terms l as i and h s i computed by the pro- 
posed procedure (solid line) and an analytical approximation 
(dashed line). The corresponding results from a long simula- 
tion run (circles) agree well with the numerical result. N = 10 5 , 
L = 10 6 , w = 0.006, p = 0.0007, r = 0.0002, . 

provide us with a first estimate cisli, cssi,i- To obtain 
a better estimate, we generate a network at the operat- 
ing point i* , Zjj, Zg S ; with third moments according to 
cisi.ij cssi.i (explained below). The endpoint of a sec- 
ond burst of simulation yields a better estimate cisi,2, 
cssi, 2- This process is repeated n- m it times to converge 
on the slow manifold (see Fig. 1). This is followed by 
Haver = ^tot — ^init further repetitions. Averaging over 
the endpoints of these repetitions we compute the consis- 
tent values c* = Y^n=n init c x,n/"avcr, for x=ISI, SSI. The 
variance of the approximation and therefore the width of 
the band in which the consistent values are located can be 
estimated by J2n=n illit ( c x ~ c x , n ) 2 / {n avcr 2 - n avcr ). 

An important step in the algorithm outlined above is the 
generation of a network with specific i, Igi, Zss> ^ssi> and 
Zisi- Starting from a given network the desired number 
of i, Zsi and hi can be reached straightforwardly by first 
flipping the state of nodes and then selectively rewiring 
links. Then we set Zssij ^isi : We chose a random link and 
additionally a random node that is not part of the link. If 
at least one of the nodes connected to the link has the same 
state as the individual node we selected, we can rewire the 
link in a way that it leaves i, Isi and hi unchanged but 
can possibly affect Zssi and Zigp If this rewiring brings us 
closer to the desired values it is accepted, otherwise a new 
node and a new link is chosen randomly. This process is 
repeated until the desired Zssi and Zisi have been reached. 

Numerical Methods. — The algorithm described 
above allows us to generate consistent states for Eq. (1) 
computationally on-demand. This enables us to use 
Eq. (1) to numerically estimate the temporal derivatives 
of the system-level state variables at a given point in state 
space. It thereby provides us with all the information that 
is needed to apply many established tools for the compu- 
tational investigation of dynamical systems. 

Let us consider a simple example. A standard task in 
the investigation of dynamical systems is the computation 
of a trajectory, starting from a given initial point xo. Per- 
haps, the most simple computational tool to perform this 
task is forward Euler integration: Given a point x n on the 
trajectory a subsequent point x, i+1 is computed by the 
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first order Taylor approximation x„ + i = x„ + r /(x)| x , 
where /(x) is the right-hand side of the equations of mo- 
tion and t is a fixed increment in time. By repeated ap- 
plication of this procedure a sequence of points is pro- 
duced that traces the desired trajectory. The approach of 
coarse-grained projective integration [22] can be used to 
apply methods like forward Euler integration to systems 
in which no closed analytic expression for the right-hand 
side of the equations of motion is available. Following this 



approach the derivatives, such as dx/dt^ 



/(x)| 



not computed analytically, but are estimated from prop- 
erly initialized bursts of microscopic simulation. 

We have performed coarse projective integration of the 
SIS model considered here. We have also used a variant 
scheme that takes some additional information from the 
moment expansion into account: In order to compute the 
temporal derivative of the state variables at a given point 
in state space x = (i,Zii,Zss)) we first use the procedure 
described in the previous section to compute appropri- 
ate Zisi and Zssi for Eq. (1). We can then close Eq. (1) 
and use it to estimate the temporal derivatives. Thus, 
in every step of this projective Euler method a new net- 
work is generated in order to compute the derivatives of 
the system-level variables. However, in order to construct 
the network and find the desired closure terms only short 
bursts of simulation are necessary; no long simulation runs 
are required. 




0.96 
0.03 



0.00, 



Time 



10k 



Fig. 3: (Color Online) Coarse projective Integration. Shown 
are time series of the infected fraction of the population, the 
density of SS-links Iss and Si-links 1st- The link numbers have 
been normalized with respect to the total number of links in 
the network. The results from a coarse projective integration 
(black) are in good agreement with a full individual-based sim- 
ulation (orange). Parameters: N = 10 5 , L — 10 6 , p = 0.06, 
r = 0.0002, p = 0.0003, n init = 2, n avcr = 1, s = 20, r = 100. 

A comparison of the coarse projective integration 
scheme with full individual-based simulation is shown in 
Fig. 3. Coarse projective integration has previously been 
proposed as a way to speed up lattice Boltzmann, ki- 
netic Monte Carlo as well as molecular dynamics simula- 
tions [22]. In the present implementation for networks, ini- 
tialization of consistent networks is time consuming, and 



thus only a small speed-up is realized (a factor of ~ 2). 
However, a significantly bigger speedup can be expected 
in larger networks. In addition there are other advan- 
tages: In the derivation of Eq. (1) the individual-based 
based stochasticity of the system is to a large extend fil- 
tered out. Furthermore, the coarse integration scheme can 
run backwards in time, which can for instance be useful in 
approximating the boundaries of basins of attraction [23] . 
But most importantly the projective Euler integration is 
a proof of principle; many more powerful methods of dy- 
namical systems theory can be applied in essentially the 
same way. 

The simple example of forward Euler integration proves 
that it is possible to close the model on-demand by us- 
ing only short bursts of microscopic simulation. This al- 
lows us for instance to use fixed-point algorithms such as 
Newton's method. These methods do not only compute 
stationary states more efficiently than integration to sta- 
tionarity; they can also be used to find dynamically un- 
stable steady states, such as saddles, that are inaccessible 
to direct simulation. 

Starting from an initial estimate x of a stationary state, 
Newton's method produces progressively better estimates 
by x n+ i = x„ — J /(x ra ), where J is the Jacobian matrix 
at x„ . As the convergence of our Newton method was not 
very sensitive to minor errors in the Jacobian we found it 
advantageous to estimate /(x) by the automated closure 
procedure and approximate the Jacobian by differentiat- 
ing the (inaccurate but faster) analytical closure. 

The Jacobian is also needed for continuation of a solu- 
tion branch and for detection of bifurcations [24] . Since 
the eigenvalues can be highly sensitive to errors in the 
Jacobian we approximate the elements of Jacobian in this 
case by finite differences of temporal derivatives computed 
numerically. Because of the sensitivity of the eigenvalues, 
averaging over many (~ 10 3 ) temporal derivative evalu- 
ations is necessary to obtain reasonably accurate eigen- 
values. A possibly better alternative is the application 
of matrix-free iterative algorithms, in which the eigenval- 
ues are determined without explicit approximation of the 
Jacobian [25]. 

Results. — By applying the Newton method for sev- 
eral parameter values, we have computed a one-parameter 
bifurcation diagram for the adaptive SIS network (Fig. 4). 
For p > 0.71 the network approaches a stationary state in 
which the prevalence of the diseases is high p ks 0.97. If p is 
lowered, a threshold is crossed (point B) at which a coarse 
limit cycle attractor emerges - most likely from a discon- 
tinuous non-equilibrium phase transition which appears 
as a fold bifurcation of limit cycles in the low dimensional 
system. We expect that a coarse saddle cycle is also born 
there. The stationary attractor and the limit cycle attrac- 
tor coexist until the stability of the stationary attractor is 
lost in a subcritical Hopf bifurcation (point A). In such bi- 
furcations an unstable limit cycle vanishes. While we were 
not able to compute this unstable oscillation directly, it is 
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Fig. 4: System-level bifurcation diagram. Branches of stable 
steady states (solid lines) and saddle points (dashed/dotted) 
have been computed by a coarse-grained Newton method. 
The white line tentatively marks a branch of saddle limit cy- 
cles. Shaded regions mark ranges of p observed during long 
individual-based simulations in the neighborhood of the at- 
tractive limit cycle (light gray) and of stable of stationary so- 
lutions (dark gray). Computation of the Jacobian eigenvalues 
reveal a subcritical Hopf bifurcation (A), two fold bifurcations 
(C,E) and a transcritical bifurcation (F). In addition there is 
a fold bifurcation of cycles (B) and a homoclinic bifurcation 
(D). Two small insets indicate the eigenvalue configuration at 
points A and C. Inset: time series on the limit cycle attractor 
at p = 0.0006. Parameters: w = 0.06, r = 0.0002, N = 10 5 , 
L = 10 6 . 

consistent with the saddle cycle mentioned above. 

At lower infection rates a stable steady state and a sad- 
dle point emerge from a fold bifurcation (point C). This 
discontinuous transition marks the edge of a narrow inter- 
val in which the limit cycle attractor coexists with the new 
stationary attractor. At slightly lower infection rates the 
limit cycle vanishes in what appears as a homoclinic bi- 
furcation (point D) . After this bifurcation the steady state 
that has emerged from C is the only attractor of the sys- 
tem. The two saddles vanish in a fold bifurcation at point 
E. Finally, at the epidemic threshold (point F), the sta- 
ble equilibrium undergoes a transcritical bifurcation with 
the trivial equilibrium (i = In = hs = 0) as it leaves the 
positive cone of the state space. 

Let us emphasize that Fig. 4 describes the dynamics on 
the emergent level: in the non-trivial stationary states the 
network is never frozen, but constantly undergoes changes 
in its topology and the states of its nodes. However, these 
changes are balanced in such a way that system-level vari- 
ables, such as i, Iss and In remain effectively stationary. 
Furthermore, note that two of the state variables (Iss 
and in) describe topological degrees of freedom. With- 
out adaptive rewiring, it is well known that the emergent 
level dynamics of the SIS model can be described by the 
single variable i and only exhibits transition E [19]. Dy- 
namics such as the self-organized cycles for SIS model are 
therefore only possible on adaptive networks. 



The dynamics of the adaptive SIS model with variable 
awareness reported here differ from the dynamics of the 
model with constant rewiring studied in [13]. Rewiring is 
in both cases a strong protective mechanism that isolates 
infected individuals and thus reduces the prevalence of the 
disease. However, in the present model this mechanism 
only plays the role if the prevalence and therefore also the 
population's awareness of the disease is high. In partic- 
ular, compared to the non-adaptive SIS model, rewiring 
does not increase the epidemic threshold (point F) in the 
present model, while a strong increase of the threshold 
was found in the model with constant rewiring. Perhaps 
the most striking difference between the model variants 
is that the awareness mechanism greatly increases the pa- 
rameter range in which oscillations can be observed. Also, 
the amplitude of the oscillations is increased considerably. 

In order to understand how variable awareness promotes 
oscillations, the underlying mechanism has to be taken 
into account. While rewiring of network connections iso- 
lates infected individuals and thus slowly decreases the 
prevalence of the disease, it also leads to the formation of a 
highly connected cluster of susceptible individuals (see [13] 
for details). If the disease manages to invade this cluster it 
can rapidly propagate causing a large outbreak. This leads 
back to a state of high prevalence completing the cycle. In 
the model with constant rewiring the invasion of the sus- 
ceptible cluster is difficult as it is still protected by rapid 
rewiring. If the amplitude of the cycle is large the number 
of infected individuals at the low point is small, which re- 
duces the chance of successful invasion further. Therefore 
large amplitude oscillations cannot persist in the model 
with constant rewiring and the parameter range in which 
oscillations are observed is small. In the model studied in 
this Letter large oscillations also result in small numbers 
of infected individuals in the low point of the cycle. How- 
ever in this case the disease is unlikely to become extinct 
as the system is still beyond the epidemic threshold and 
the protection offered by rewiring decreases proportional 
to the number of infected. Since the disease free state is 
thus dynamically unstable reinfection of the susceptible 
cluster is bound to occur. 

Summary and Conclusions. — In this Letter we 
have proposed a computational approach to the investi- 
gation of emergent properties of adaptive networks. By 
applying the methods of coarse-grained, 'equation-free' 
modeling and analysis we have shown that existing nu- 
merical tools of dynamical systems theory can be used to 
investigate networks directly on an emergent level. The 
approach enabled us to obtain information that is not di- 
rectly accessible by simulation alone, while avoiding the 
often prohibitively difficult derivation of closed system- 
level equations of motion. In particular we have avoided 
the strong homogeneity assumption that is inherent in pre- 
vious analytical moment closure approximations. 

Our approach depends critically on the assumption that 
bona fide emergent-level variables exist and can be iden- 
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tificd; else closure (numerical or otherwise) at the desired 
level is impossible. The moment-based description used 
here is applicable to networks in which the number of 
different states of the nodes is small and the topological 
change can be expressed in terms of the network moments. 
The success of this description depends critically on the 
time scale separation between the state variables and the 
higher moments of the network. In our model the time 
scale separation is finite and decreases as the rewiring rate 
is increased. At very high rewiring rates the computation 
of closure terms ceases to converge. 

In other systems it may be advantageous to use different 
sets of system-level variables, such as heuristic variables, 
based on the researcher's experience or variables identi- 
fied by automated data- mining algorithms [26]. Although 
dynamical equations, such as Eq. (1) are in this case not 
available, one can still estimate temporal derivatives of 
the chosen set of variables directly from short bursts of 
properly initialized simulation. 

Here we only applied two simple numerical tools: For- 
ward Euler integration and the Newton method. These 
tools revealed significant information on the coarse bifur- 
cation structure of the system. Demonstrating the ap- 
plicability of such tools to a network problem constitutes 
a proof of principle. More sophisticated integrators and 
fixed point algorithms, continuation of solution branches, 
automated bifurcation detection, computation of normal 
form coefficients etc. can in principle be analogously ap- 
plied. These tools are available in the form of free, well- 
tested software packages that are commonly used for the 
investigation of systems of ordinary differential equations. 
By writing numerical wrappers that operate along the 
lines discussed above, this arsenal of existing, highly effi- 
cient methods can be brought to bear on the investigation 
of adaptive networks. 

In the present Letter we have used the proposed proce- 
dure to investigate an epidemiological SIS model. In our 
model individuals can protect themselves by altering their 
topological neighborhood at a rate determined by their 
awareness to the disease. Our analysis shows that in a 
considerable parameter range the prevalence of the dis- 
ease and the topology of the network exhibits oscillations 
of large amplitude. Let us emphasize that the observed 
oscillations are not caused by the variable level of aware- 
ness, since in our model awareness changes without time 
lag. Instead variable awareness levels extend the parame- 
ter range in which oscillations can be observed; decreasing 
awareness levels at low prevalence stabilize cycles which 
would otherwise drive the disease to extinction. The sim- 
ple mechanism studied in our conceptual model may ex- 
plain oscillatory dynamics observed in certain real world 
diseases. 
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